#################################################################
##DO SELF-REPORTING REGIMES MATTER? EVIDENCE FROM THE CAT########
## Cosette D. Creamer and Beth A. Simmons #######################
## February 2019 ###############################################
### Replication Code for Appendix B ##########
### ########################### ################################
#################################################################
#################################################################
#################################################################
## Clear environment
rm(list=ls())
sink("CAT_AppendixB.txt")
### Libraries
library(data.table)
library(MASS)
library(car)
library(SparseM)
library(stargazer)
library(Zelig)

##############################################
### Data loading and pre-processing  ##########
###############################################
## Load CAT data (full)
catfull <- read.csv("CATReporting_Final.csv", header=TRUE)

## Subset data to only include observations for CAT ratifiers, for observations > t+1 (t = year of ratification)
cat <- subset(catfull, catfull$CATob==1)
cat$CATregiondensityl1 <- 100*cat$CATregiondensityl1
cat$HRTpercentl1 <- 100*cat$hrtpropl1
cat$yr <-as.factor(cat$year)
cat$lngdpl1 <- log(cat$gdpl1)
cat$lnpop <- log(cat$poptot)

## Remove observations with missingness for covariates used 

cat <- subset(cat, !is.na(cat$lngdpl1) & !is.na(cat$lnpop) & !is.na(cat$tortl1) & !is.na(cat$farissl2) &
                    !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
                  & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) &
                    !is.na(cat$regionall1) & !is.na(cat$transdeml1) & 
                    !is.na(cat$transition2l1) &!is.na(cat$HRTpercentl1) &!is.na(cat$transparencyindex_l1))


########################################################################
# Denominator model (YEAR + treatment(review) history variables + time-varying covariates)
denom.fit <-  glm(CATreview ~ lngdpl1 + lnpop 
                  +nhril1 +transparencyindex_l1 +polityl1+ nevdem + transdeml1+transition2l1+ 
                    farissl2+ tortl1 +HRTpercentl1 + CATart22l1
                  +   CATregiondensityl1+regionall1+
                    CATnumprevrevs+CATrevprevyr + CATyrssincereview + 
                    CATlintt+farissl2*CATlintt, 
                  family="binomial", data=cat)
denom.p <- predict(denom.fit, type = "response")

#################################################################
### Table B1: Estimated Probability of undergoing CmAT Review  ##
# Numerator model (YEAR + treatment(review) history variables) ##
#################################################################
stargazer(denom.fit)
#################################################################
#################################################################
numer.fit <-  glm(CATreview ~ 
                    CATrevprevyr + CATyrssincereview + CATnumprevrevs + nevdem+ CATlintt, 
                  family="binomial", data=cat)
numer.p <- predict(numer.fit, type = "response")
## Calculate weights


d.pr.tr <- ifelse(cat$CATreview == 1, denom.p, 1-denom.p)
n.pr.tr <- ifelse(cat$CATreview == 1, numer.p, 1-numer.p)
cat$w <- n.pr.tr/d.pr.tr
summary(cat$w)
sd(cat$w)
dt <- data.table(cat)
dt[ , swrev:=NA_real_]
dt[ intersect(order(year), which(!is.na(w))), swrev := cumprod(w), by=cowcode]
cat <- as.data.frame(dt)
summary(cat$swrev)
sd(cat$swrev)



############################################################
#### Outcome MSM with CIRI Torture Score (t+2)  #######
############## Table B2   #####################
############################################################
cat$tortf2 <- factor(cat$tortf2)
cat2 <- subset(cat, !is.na(cat$tortf2))
msm5 <- polr(tortf2 ~ CATreview + CATnumprevrevs+  CATrevprevyr+ CATyrssincereview+ 
       farissf1+CATlintt+farissf1*CATlintt ,method= 'probit',  weights = swrev, 
     data = cat2)

summary(msm5)
stargazer(msm5)
#############################################################################
#########Figure B2###########################################################
### Note: this uses output from analysis in "CAT_CIRI MSM Analysis.R "########
##############################################################################
model3<-read.csv("CAT_CIRI__model_3_alpha_0_iterations_20000_outcomelead_2.csv", header=TRUE)


m3_coeff<-rbind(model3$CATreview,model3$CATrevprevyr,model3$CATyrssincereview,model3$CATnumprevrevs,model3$farissf1,model3$CATlintt,model3$farissf1.CATlintt)

m3_95ci_upper <- apply(model3, 2, function(x) quantile(x, .975))
m3_95ci_lower <- apply(model3, 2, function(x) quantile(x, .025))
m3_mean<- apply(model3, 2, function(x) quantile(x, .50))
out3mean<-(m3_mean)
out3low<-(m3_95ci_lower)
out3upper<-(m3_95ci_upper)
g3<-rbind(out3mean,out3low,out3upper)
ci3<-t(g3)
stargazer(ci3)

xc.nrrev<-cbind(0, 0, 2.63, c(0:5), .45, 13.2, 5.94)
head(xc.nrrev)

ystar.nrrev<-(xc.nrrev)%*%(m3_coeff)
ystar.nrrev<-t(ystar.nrrev)
dim(ystar.nrrev)
m3_zeta<-cbind(model3$X0.1, model3$X1.2)
m3_zeta<-as.data.frame(m3_zeta)
pr1<- pnorm(m3_zeta$V1,mean=ystar.nrrev,sd=1)
pr2<- pnorm(m3_zeta$V2,mean=ystar.nrrev,sd=1)-pnorm(m3_zeta$V1,mean=ystar.nrrev,sd=1)
pr3<-1-pnorm(m3_zeta$V2,mean=ystar.nrrev,sd=1)
pr1mean<-apply(pr1,2,median)
pr2mean<-apply(pr2,2,median)
pr3mean<-apply(pr3,2,median)


postscript("FigureA8.eps", width=5.5, height=4.5, horizontal = FALSE, 
           onefile = FALSE)


par(cex.axis=1, cex.lab=0.9, cex.main=1, cex.sub=1)
plot(c(0:5),pr1mean,ylim=c(0,1),xlim=c(0,5),lwd=2,type="l",lty=3,
     xlab="Number of Previous Reviews ", ylab="Probability of observing a given level of Torture")

lines(c(0:5),pr2mean,ylim=c(0,1),lwd=2,lty=2,type="l")
lines(c(0:5),pr3mean,ylim=c(0,1),lwd=2,lty=4,type="l")
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(3,0,0),lwd=2, bty="n",cex=.7)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,2,0),lwd=2, bty="n",cex=.7)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,0,4),lwd=2, bty="n",cex=.7)

p1.low<- apply(pr1,2,quantile,prob=.025)
p1.high<- apply(pr1,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p1.low,  
         y1=p1.high,lwd=1) 
p2.low<- apply(pr2,2,quantile,prob=.025)
p2.high<- apply(pr2,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p2.low,  
         y1=p2.high,lwd=1) 
p3.low<- apply(pr3,2,quantile,prob=.025)
p3.high<- apply(pr3,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p3.low,  
         y1=p3.high,lwd=1)
dev.off()

sink(NULL)
